MODULE ModCoolCII

!!$ This module contains the subroutine for calculating the cooling due to CII.
!!$ Programming conventions are from Nyhoff and Leestma, _Fortran 90 for Scientists and Engineers_, 1997
!!$
!!$  STH: Samuel Harrold

!!$ Declare modules to USE
  USE ModParms
  USE ModFuncs
  USE ModAbund

!!$ Do not assume undeclared variables
  IMPLICIT NONE



CONTAINS



!!$ CalcCoolCII
!!$ Calculate the cooling due to the magnetic dipole fine-structure forbidden transition [CII] 2P(3/2) -> 2P(1/2) at 158 microns for given physical conditions. The cooling function follows the 2-level treatment described in Tielens and Hollenbach 1985 (TH85).
!!$ Accepts: particle number density, temperature, mean intensity of radiation field with respect to cm^-1, column density in line of sight to the star
!!$ gg frequency index nu for 158 microns is 13 because
!!$ lambda = 10^3 10^(-4n/64) in microns for n = 1,2,...64
!!$ TODO: calclate by arrays rather than by element

  SUBROUTINE CalcCoolCII(ngas, Tempgas, Jw, Ntostar, CoolCII)

!!$ Declare variables
!!$ Specify input and output
!!$ ngas    = density of gas in particles cm^-3
!!$ Tempgas = temperature of gas in K
!!$ Jw      = mean intensity of radiation field in erg/(s cm^2 sr cm^-1)
!!$ w for omega as wavenumber
!!$ Ntostar = column density of gas in line of sight to star in particles cm^-2
!!$ CoolCII = cooling due to emission from CII in units acceptable by pisco 
!!$ in erg s^-1 (gas particle)^-1
    REAL, INTENT(IN)  :: ngas , Tempgas, Jw, Ntostar
    REAL, INTENT(OUT) :: CoolCII

!!$ Specify internal variables
!!$ Data for [CII] at 158 microns from NIST Atomic Spectra Database
!!$ Aul  = radiative de-excitation coefficient in s^-1
!!$ glwr    = statistical weight g lower
!!$ gupr    = statistical weight g upper
!!$ nu0     = frequency of CII transition at 157.6790 microns in Hz
    REAL :: &
         Aul  = 2.29E-6, &
         glwr = 2., &
         gupr = 4., &
         nu0  = c / (157.6790 * 1.E-4)

!!$ Variables for physical conditions
!!$ Notation follows that of Tielens and Hollenbach 1985 (TH85).
!!$ TUarraye = array of coordinate paris of gas temperature (T) and CII-e de-excitation collision strength (U for Upsilon). Approximated from Blum and Pradhan 1992 and Osterbrock book 2006, p51. Note: only valid for 1 <= Tempgas <= 5000
!!$ Jnu = mean intensity of radiation field wrt Hz, units erg/(s cm^2 sr Hz)
!!$ Abund(CII/e/H/H2) = abundances of CII/e/H/H2 relative to all gas particles determined by column density to the star
!!$ n(CII/e/H/H2)    = number density of CII (all levels)/e/H/H2
!!$ CollisionStre = CII-e excitation collision strength
!!$ Besc     = escape probability
!!$ Pnu      = background radiation field wrt Hz
!!$ Qnu      = population mode of background radiation field wrt Hz
!!$ gamma(ul/lu)(e/H/H2) = collisional excitation/de-excitation. Collision strength formulafrom Osterbrock book 2006, p51. H from Launay, Roeff 1977b via TH85. H2 from Flower, Launay 1977 via TH85.
!!$ C(ul/lu)   = composite collisional excitation/de-excitation coefficient
!!$ R(ul/lu)   = rate of transition from upper to lower/lower to upper excitation state
!!$ n(upr/lwr) = population density of upper/lower state
!!$ Snu = source function
!!$ L = cooling from TH85, Eqn B1 in units of erg s^-1 cm^-3
    REAL, DIMENSION(5, 2) :: TUarraye
    REAL :: &
         Jnu, &
         AbundCII, Abunde, AbundH, AbundH2, &
         nCII, ne, nH, nH2, &
         CollisionStre, &
         Besc, &
         Pnu, Qnu, & 
         gammaule,  gammalue, &
         gammaulH,  gammaluH, &
         gammaulH2, gammaluH2, &
         Cul, Clu, &
         Rul, Rlu, &
         nupr, nlwr, &
         Snu, &
         L

!!$ Convert mean intensity of radiation field
!!$ Jnu = mean intensity of radiation field converted to erg/(s cm^2 sr Hz)
!!$ using nu = c*w => dnu = c*dw and Jnu is defined as a derivative
!!$ wrt nu, Jnu = dJ/dnu = dJ/dw * dw/dnu = Jw/c
    Jnu = Jw / c

!!$ Estimate number densities of species
!!$ number density of species = abundance * number density of gas
    CALL CalcAbund(Ntostar = Ntostar, &
         AbundCII = AbundCII, Abunde = Abunde, &
         AbundH = AbundH, AbundH2 = AbundH2)
    nCII = ngas * AbundCII
    ne   = ngas * Abunde
    nH   = ngas * AbundH
    nH2  = ngas * AbundH2

!!$ Estimate collision strength for CII-e collisional excitations.
!!$ Source: Blum, Pradhan 1992; Osterbrock, Ferland 2006.
!!$ Note: Only valid for 1 <= Tempgas <= 5000
!!$ Array of temperature and velocity-averaged collision strengths TUarraye used to interpolate collision strength for a given temperature CollisionStre.
!!$ InterpLinear interpolates a linear function given two coordinate pairs:
!!$ y = (x - xi) * (yf - yi)/(xf - xi) + yi
    TUarraye = RESHAPE( (/ &
         1000., 1.5763, &
         2000., 1.6361, &
         3000., 1.7157, &
         4000., 1.8034, &
         5000., 1.8874 /), (/5, 2/), &
         ORDER = (/2, 1/) )
    IF (Tempgas < TUarraye(2, 1)) THEN
       CollisionStre = InterpLinear(x = Tempgas, &
            xi = TUarraye(1, 1), yi = TUarraye(1, 2), &
            xf = TUarraye(2, 1), yf = TUarraye(2, 2))
    ELSE IF (TUarraye(2, 1) <= Tempgas .AND. Tempgas < TUarraye(3, 1)) THEN
       CollisionStre = InterpLinear(x = Tempgas, &
            xi = TUarraye(2, 1), yi = TUarraye(2, 2), &
            xf = TUarraye(3, 1), yf = TUarraye(3, 2))
    ELSE IF (TUarraye(3, 1) <= Tempgas .AND. Tempgas < TUarraye(4, 1)) THEN
       CollisionStre = InterpLinear(x = Tempgas, &
            xi = TUarraye(3, 1), yi = TUarraye(3, 2), &
            xf = TUarraye(4, 1), yf = TUarraye(4, 2))
    ELSE IF (TUarraye(4, 1) <= Tempgas) THEN
       CollisionStre = InterpLinear(x = Tempgas, &
            xi = TUarraye(4, 1), yi = TUarraye(4, 2), &
            xf = TUarraye(5, 1), yf = TUarraye(5, 2))
    END IF

!!$ Calculate escape probability, Besc
!!$ see de Jong, Dalgarno, and Boland 1980 for plane parallel approximation
!!$ Assume escape probability is 0.5 from face of disk
!!$ TODO: make tauIR an input
    Besc = 0.5

!!$ Next calculations follow Tielens and Hollenbach 1985 (TH85)
!!$

!!$ Calculate background radiation field for each transition wrt Hz in erg / (s cm^2 sr Hz) and the  mode of radiation field for each transition (unitless)
    Pnu = Jnu
    Qnu = (c**2 / (2*h*nu0**3)) * Pnu

!!$ Calculate collisional coefficient for excitation and deexcitation for each transition and species. Listed from most to least important
!!$ Collision coefficient gamma e from Osterbrock, Ferland 2006; Blum, Pradhan 1992
    gammaule  = 8.629E-6 * (CollisionStre/gupr) * Tempgas**(-0.5)
    gammalue  = (gupr/glwr) * gammaule  * EXP(-h*nu0 / (k*Tempgas))

!!$ gamma H from Launay and Roeff 1977b via Tielens, Hollenbach 1985
    gammaulH  = 5.8E-10 * Tempgas**0.02
    gammaluH  = (gupr/glwr) * gammaulH  * EXP(-h*nu0 / (k*Tempgas))

!!$ gamma H2 from Flower and Launay 1977 via Tielens, Hollenbach 1985
    gammaulH2 = 3.1E-10 * Tempgas**0.1
    gammaluH2 = (gupr/glwr) * gammaulH2 * EXP(-h*nu0 / (k*Tempgas))

!!$ Calculate collisional excitation/de-excitation rates, units transition/(s cm^3). Terms are listed from most to least important.
    Cul  = gammaule*ne + gammaulH*nH + gammaulH2*nH2
    Clu  = gammalue*ne + gammaluH*nH + gammaluH2*nH2

!!$ Calculate excitation/de-excitation rates, units transition/(s cm^3)
    Rul  = Aul*Besc*(1. + Qnu) + Cul
    Rlu  = (gupr/glwr)*Aul*Besc*Qnu + Clu

!!$ Calculate level populations populations, nupr, nlwr
    nupr = nCII * (1. + (Rul/Rlu))**(-1)
    nlwr = nCII - nupr

!!$ Calculate source function, units erg/(s cm^2 sr Hz)
    Snu = (2.*h*nu0**3/c**2) * (( (gupr*nlwr)/(glwr*nupr) ) - 1.)**(-1)

!!$ Calculate the cooling due to CII, CoolCII
!!$ Pisco needs cooling in erg s^-1 (gas particle)^-1
!!$ From TH85, Eqn B1, Lambda is in erg s^-1 cm^-3
!!$ Cooling is thus Lambda/ngas (ngas = particle density of gas input)
    L = nupr*Aul*h*nu0*Besc*(1. - (Pnu/Snu))
    CoolCII = L/ngas

  END SUBROUTINE CalcCoolCII



END MODULE ModCoolCII
